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ABSTRACT 

Numerical  bifurcation  techniques  were  developed  for  studying 
the  multiplicity,  stability,  and  oscillatory  dynamics  of  the  non- 
adiabatic  tubular  reactor  with  a  single  A  B  reaction.  The 
techniques  illustrate  the  existence  of  one,  three,  five,  or  seven 
steady  states  and  bifurcating  periodic  solutions.  We  present 
numerical  procedures  for  computing  the  Hopf  bifurcation  formulas 
which  can  determine  the  stability  and  location  of  the  oscillation 
without  integrating  the  parabolic  partial  differential  equations. 
The  combination  of  our  Hopf  techniques  with  steady  state  bifurca¬ 
tion  methods  enables  us  to  determine  all  possible  steady  and  stable 
oscillatory  solutions  exhibited  by  distributed  parameter  models 
such  as  the  tubular  reactor. 
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SIGNIFICANCE  AND  EXPLANATION 


Numerical  methods  are  developed  for  the  investigation  of  Hopf  bifurcation 
(bifurcation  to  periodic  solutions)  in  mathematical  models  which  consist  of 
parabolic  partial  differential  equations  whoso  time  independent  solutions 
are  defined  by  systems  of  two  point  boundary  value  problems.  The  combination 
of  our  Hopf  techniques  with  the  steady  state  bifurcation  methods  of  Keller 
enables  us  to  determine  all  possible  steady  state  and  periodic  solutions 
exhibited  by  these  distributed  parameter  models.  The  utility  of  our  numerical 
procedure  lies  in  their  generality  and  potential  applicability.  They  can  be 
used  to  study  the  multiplicity,  stability,  and  oscillatory  phenomena  exhibited 
by  reaction-diffusion  systems  in  combustion,  chemical  reactor  theory  and  mathe¬ 
matical  biology. 

We  apply  these  techniques  to  a  model  of  the  tubular  reactor  with  an 
exothermic  reaction.  The  reactor  is  found  to  exhibit  broad  regions  of  multiple 
steady  states  and  periodic  solutions  which  can  be  conveniently  presented  on 
response  curves  to  illustrate  the  dynamic  capabilities  of  this  reactor.  We 
discuss  the  effects  of  the  solutions  on  reactor  operation,  and  the  effects  of 
sharp  transitions  or  jumps  between  the  steady  states  and  oscillatory  states  on 
reactor  dynamics.  Many  of  the  above  results  have  not  been  reported  in  the 
earlier  studies  of  the  tubular  reactor  or  in  the  more  exhaustive  studies  of 
the  continuously  stirred  tank  reactor. 


T'ne  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  li*  s  wit.li  MRC,  and  not  with  the  authors  of  this  report. 


NUMERICAL  HOPF  BIFURCATION  TECHNIQUES 
AND  THE  DYNAMICS  OF  THE  TUBULAR  REACTOR  MODEL 


Robert  F.  Heinemann  and  Aubrey  B.  Poore 
INTRODUCTION 

The  number  of  theoretical  and  experimental  investigations  of  multiplicity, 
stability,  and  sensitivity  of  reaction-diffusion  systems  over  the  past  two 
decades  has  been  enormous.  Much  of  this  research  has  centered  around  heat 
and  mass  transfer  coupled  with  chemical  reaction  and  has  generally  been 
motivated  by  the  design  and  operation  of  industrial  chemical  reactors.  The 
most  widely  studied  reactor  type  has  been  the  continuously  stirred  tank  reactor 
(CSTR)  since  the  relevant  mathematics  are  comparatively  tractable,  and  studies 
of  the  CSTR  provide  a  foundation  for  the  understanding  of  more  complex  reactors. 
Poore  [1]  and  Uppal  et  al.  [2,3]  have  presented  an  important  and  thorough  con¬ 
sideration  of  the  single,  first-order  exothermic  reaction  in  the  CSTR.  They 
concisely  classified  all  possible  patterns  of  multiple  steady  and  periodic 
states  and  pieced  together  a  myriad  of  results  presented  in  earlier  papers. 

The  Hopf  bifurcation  techniques  presented  by  Poore  and  Uppal  et  al .  were  com¬ 
pletely  general  for  time-dependent  ordinary  differential  equations  and  have  been 
applied  to  other  systems  exhibiting  oscillatory  phenomena  [4] .  Such  systems 
included  those  whose  mathematical  models  consist  of  parabolic  partial  differ¬ 
ential  equations  which  can  be  judiciously  simplified  to  ordinary  differential 
equations  [5] .  These  simplifications  have  been  advantageous  since  Hopf  bifur¬ 
cation  techniques  were  not  available  for  distributed  parameter  models.  The 
major  objective  of  the  present  work  is  to  fill  this  void  by  presenting  the 
Hopf  formulas  for  these  complex  systems. 

Our  numerical  techniques  are  based  entirely  on  the  time-independent 
problem  and  yield  the  relevent  bifurcation  information  such  as  the  direction 
and  stability  of  the  oscillatory  solution.  Furthermore,  the  periodic  orbit 
Sponsored  by  the  United  States  Army  under  Contract  No.  DAAG29-80-C-0041 . 
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can  be  computed  without  time  integration  from  an  asymptotic  formula  for  the 
orbit.  These  methods  are  quite  general  and  can  be  applied  to  systems  of 
parabolic  partial  differential  equations  whose  time-independent  solutions  are 
defined  by  a  two-point  boundary  value  problem. 

We  illustrate  the  utility  of  our  numerical  techniques  by  applying  them  to 
the  model  of  the  nonadiabatic  tubular  reactor.  Because  of  its  fundamental 
importance,  this  reactor  has  been  studied  with  only  slightly  less  fervor  than 
the  CSTK.  The  formidable  body  of  literature  concerned  with  the  problem  includes 
the  extensive  work  of  Amundson  et  al.  [5  -  15],  Hlavecek  et  al .  [16  -  19],  and 

McGowin  and  Perlmutter  [20].  Many  of  the  results  presented  in  these  works  are 
summarized  and  explained  in  two  excellent  reviews  by  Schmitz  [21]  and  Varma  and 
Aris  [22],  These  earlier  papers  have  firmly  established  the  existence  of  one, 
three,  or  five  steady  states  and  also  have  illustrated  sustained  oscillations. 
However,  a  complete  understanding  of  the  oscillatory  dynamics  is  still  largely 
an  un t one h e d  problem. 

Two  recent  papers  concerning  the  multiplicity  of  the  nonadiabatic  tubular 
reactor  should  be  mentioned  at  this  point.  Kapilla  and  Poore  [23]  have  completely 
c! ct:  sir ied  the  structure  of  the  multiple  steady  states  for  all  possible  para- 
Hcfa's  using  large  activation  energy  asymptotics  and  have  established  the  exist¬ 
ence  of  two  additional  solutions  not  seen  in  previous  works*  We  [24]  have  con¬ 
firmed  and  extended  this  classification  by  applying  the  steady  state  bifurcation 
techniques.;  of  Keller  to  the  problem  [25],  The  effect  of  a  wide  range  of  finite 
activation  energies  has  been  examined  and  has  shown  the  existence  of  one,  three, 
five  and  seven  solutions.  Following  presentation  of  th^  mathematical  model  of  the 
re-actor  and  our  numerical  procedures,  we  illustrate  Hopf  bifurcation  results  for 


multiplicity  patterns  exhibiting  from  or.o  to  seven  steady  states. 


MATHEMATICAL  MODEL 


The  equations  describing  the  conservation  of  reactant  A  and  energy  for 
the  nonadiabatic  tubular  reactor  with  axial  mixing  appear  below  in  dimension¬ 


less  form: 


i£=-i_9_Z_iX_  Dv  e>_Y/y 

3t  Pe  .  2  3s  y 
m  3s 

3  0  1  3^6  30  „  t  .  x  y-y/0 

T-  =  - -  - g  -  T - £  +  B  Dy  e1  r/ 

3x  Pe  .  2  3s  0 

h  3s 

The  boundary  and  initial  conditions  are: 


Pe  (y-1)  at  s  =  0  ,  t  >  0 
m 


t —  -  Pe  (0-1)  at  s  =  0  ,  t  >  0 
3s  h 


3y  _  30 
3s  3s 

y  =  y. 


at  s  =  1  ,  t  >  0 


e  =  e.  at  t  =  o  . 
m 


In  writing  these  equations,  we  have  defined  the  following  dimensionless 


quantities  , 


y  =  c/cr 


S  =  X/L 


Pe  =  vL/D 
m  e 


B  - AH  Wo 

D  =  Ae  ^L/v 


0  =  T/TQ 
r  =  tv/L 


Pe,  =  pC  vL/k 
h  p  e 

B  =  UPL/avpC 


y  =  E/RTQ  . 


The  above  model  describes  an  exothermic  A  ►  B  reaction  occurring  in  a  homoge¬ 
neous  tube  under  the  assumptions  that  the  velocity  profile  is  flat  with  constant 
velocity  v  ;  the  variables  y  and  °  depend  only  on  one  space  dimension  and 
time;  the  diffusion  of  reactant  A  is  governed  by  Pick's  Law  with  an  effective 
diffusivity,  D  ;  heat  conduction  is  described  by  Fourier’s  Law  with  an  effective 


thermal  conductivity,  ;  the  heat  loss  at  any  point  is  proportional  to 

(9  -  6q)  j  and  the  reaction  rate  is  describable  by  an  Arrhenius  expression. 

The  dimensional  predecessors  of  the  above  equations  and  the  applica¬ 
bility  of  this  formulation  are  discussed  in  detail  in  the  two  review  article 
and  in  the  earlier  reactor  papers. 
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NUMERICAL  METHODS 


We  now  present  the  numerical  techniques  for  generating  the  steady  state 
response  curves  and  the  Hopf  bifurcation  information  for  systems  whose  mathe¬ 
matical  description  is  given  by  a  distributed  parameter  model.  The  coupling 
of  these  methods  enables  us  to  systematically  locate  all  the  possible  steady 
and  periodic  states  exhibited  by  these  systems.  We  begin  with  a  discussion 
of  the  general  forms  of  the  equations  that  we  are  capable  of  treating. 

Many  of  the  distributed  parameter  models  in  chemical  reactor  theory  can 

be  written  as  a  system  of  parabolic  partial  differential  equations  of  the  form 

=  A (u )  - g-  +  f  (u,  u  #  u  ) 

3t  .2  x 

3x 

G  (u  (a)  ,  u  (a)  ,u)  =  0  (6) 

i  x 

G  (u  (b)  ,  u  (b)  ,  jj)  =  0 
r  x 

where  u  may  represent,  for  example,  the  concentration  and  temperature  in  the 
reactor.  G,,  and  G  denote  the  left  and  right  boundary  conditions  (possibly 
nonlinear)  and  p  is  the  bifurcation  parameter.  The  function  f  may  depend 
nonlinearly  on  u,  u^,  and  p  and  may  contain  the  space  variable  x  explic¬ 
itly  but  not  the  time  variable  t  .  For  brevity  we  write  this  system  as 
=  F ( u ,  p )  ,  G(u,p)  =  0  .  (7) 

The  corresponding  two-point  boundary  value  problem  is  denoted  by 

F(v,p)  =  0  ,  G (v , u )  =  0  .  (B) 

Steady  State  Numerical  Techniques 

The  steady  state  problem  (8)  was  solved  by  combining  Keller's  modification 
of  the  Euler-Newton  continuation  method  [25]  with  de  Boor  and  Wiess'  spline 
collocation  code  (26]  and  a  fourth-order  finite  difference  scheme  due  to  Stpplei 
[27].  Both  the  collocation  and  finite  difference  methods  performed  quite  well 
with  the  collocation  technique  capable  of  higher  accuracy. 

The  objective  of  the  steady  state  methods  is  to  compute  the  solution  v 
of  equation  (8)  as  the  parameter  u  assumes  all  its  possible  values.  The 


Euler-Newton  continuation  method  can  be  used  to  solve  this  problem  with  Euler's 

method  (9)  serving  as  a  predictor  for  Newton's  method  (10) 

v°U.  +  6u>  =  (9) 

a  u 

V^  +  ^(Vi  +  i\.)  =  v^(u+cjj)  -  F  (v*  (u+5u)  ,  U+Su)  F  (v*  (vi+5v- )  /  u+^u)  -  (10) 

However,  the  technique  fails  near  transition  between  steady  states  since  the 
Jacobian  matrix,  F  ,  cannot  be  inverted  near  singularities  such  as  limit  or  bi¬ 
furcation  points. 

Keller  has  modified  the  above  method  by  imposing  an  additional  normalization 
on  the  solution  which  enables  entire  solution  branches  to  be  traced,  skipping 
over  any  singular  points.  The  imposition  of  this  constraint  allows  the  specifi¬ 
cation  of  a  new  parameter,  s  ,  which  replaces  u  as  the  continuation  parameter 
in  the  Euler-Newton  technique.  The  reparameterized  problem  becomes 

P(x,s)  =  0  (11) 

where 

x  (s )  =  (v(s) , l  vs) )  (12) 

and 

j  F ( v (s )  ,  ^  (s)  )  j 
P (x(s) ,s)  =  )  I 

{  N  (v  (s)  ,  (s)  ,s)  j 

It  is  convenient  to  choose  the  normalization  N(v,u,s)  so  that  s  approximates 
the  arc-length  of  the  solution  branch  for  some  parameter  •  (0,1) 

N(v,,,s)  =*  «l!  v  (s)  -  v(sQ)ir  +  {l-.0|u(s)  -  t(s0)|2  -  (s  -  sQ)2.  (13) 

When  the  Euler-Newton  technique  is  applied  to  (11),  computational  diffi¬ 
culties  near  singularities  are  eliminated  since  the  Jacobian  matrix  remains  non¬ 
singular  near  such  points. 

These  techniques  and  an  algorithm  for  their  implementation  are  presented 
in  detail  by  Keller  125]  and  have  previously  been  applied  to  laminar  flame  problems 
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[28-30]  and  catalysts  problems  [31]  which  exhibit  multiple  steady  states.  There¬ 
fore,  we  forego  a  discussion  of  the  exact  computational  procedure. 

Hopf  Bifurcation  Formalism 

We  use  the  term  formalism  here  since  the  presentation  will  be  stripped  of 
the  technical  mathematical  assumptions.  A  proper  mathematical  framework  can  be 

found  in  the  works  of  Crandall  and  Rabinowitz  [32]  and  loos  and  Joseph  [33]. 

Our  work  follows  closely  the  presentation  in  the  latter  paper,  though  modified 
somewhat  to  account  for  the  nonzero  steady  state  problem  and  the  form  of  the 

model  equations  (7) .  Since  this  bifurcation  theory  is  most  effectively  used 
in  a  study  of  the  dynamics  associated  with  an  exchange  of  stability,  we  begin 

with  a  brief  discussion  of  steady  state  stability. 

The  stability  of  time-independent  solutions  can  in  principle  be  resolved  by 
examining  the  eigenvalues  of  the  linearized  boundary  value  problem.  If  the 
eigenvalues  all  have  negative  real  parts,  the  steady  state  is  stable;  whereas,  if 
an  eigenvalue  has  a  positive  real  part,  the  solution  is  unstable.  Exceptions  to 
this  principle  occur  when  the  linearized  problem  has  a  zero  eigenvalue  or  a  pair 
of  complex  conjugate,  purely  imaginary  eigenvalues.  In  the  current  reactor 
problem,  a  zero  eigenvalue  gives  rise  to  a  limit  point  bifurcation  (a  point  of 
vertical  tangency)  on  the  response  curves.  The  bifurcation  of  a  periodic 
solution  (Hopf  bifurcation)  occurs  when  a  pair  of  complex  conjugate  eigenvalues 
j(w)  and  o(u)  become  purely  imaginary.  We  assume  that  this  crossing  of  the 
imaginary  axis  occurs  at  so  that  o ( U Q with  w  positive.  It  is 

also  assumed  that  Re  o  (  i»  )  f  0  where  o  =  do/d;i.  This  ensures  a  strict 

is  0 

crossing  of  the  axis  and  is  nearly  always  satisfied  in  these  problems. 


When  the  Hopf  point  is  accompanied  by  an  exchange  of  stability,  different 


types  of  periodic  phenomena  may  be  encountered.  If  the  periodic  orbit  is 
stable,  a  small  amplitude  oscillation  is  observed  near  the  Hopf  point,  but  if 
the  orbit  is  unstable,  the  solution  will  jump  to  either  a  large  amplitude, 
stable  oscillation  or  to  another  stable  steady  state.  Examination  of  the  orbit 
stability  near  the  Hopf  points  suggests  a  systematic  procedure  for  locating  the 
stable  oscillations  in  the  reactor. 

To  present  the  Hopf  bifurcation  formulas,  we  linearize  the  boundary  value 
problem  and  write  it  as 

Lw=0,  G(vU,^|w)=0  (14) 

U  v  1 

where 

L  W  =  F  (vU,p I W)  =  ~  F (VU  +  6w,u) | 

jj  V  d  u  o- ■  0 

and 

G  (vki  ,p  |w)  =  G  (vU  +  6w,u)  |  . 

V  do  0SU 

The  essential  requirements  for  Hopf  bifurcation  without  their  technical  assump¬ 
tions  may  be  summarized  as  follows.  Assume  lico^  are  simple  eigenvalues  of  L  , 
that  niw^  is  not  an  eigenvalue  for  n  =  0,2,3,...,  and  that  the  real  part  of 
e  (Uq)  is  nonzero.  Then  one  can  construct  a  bifurcating  periodic  solution  of  (7) 
with  frequency  wj(e)  via  a  perturbational  expansion  which  can  be  shown  to  take 


the  form  [34] 

u (x, t) 


+ 


ic 

=  v  +  2  e  Re  {  Cq  e  } 

2  P0 

~  {wL  +  — u  +  2  Re{w2e2iS}}  +  0(c2) 

2 

+  \  -2  +  0  (£2) 


w  (  t) 


+  0(c2) 


(15) 

(16) 
(17) 


t  =  t.'  ( r. )  s 


(18) 
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c  is  an  auxiliary  parameter  representing  the  amplitude  of  the  orbit  close-  to 
the  Hopf  point.  The  vector  function  £  is  the  eigenvector  corresponding  to 
the  eigenvalue  +  i^  ;  and  W2  are  s°luti°ns  certain  linear  nonhomo- 

geneous  boundary  value  problems  discussed  in  the  next  section.  We  note  that 
the  sign  of  \x  yields  the  sign  of  u  -  for  £  sufficiently  small  and 

therefore  determines  the  direction  of  bifurcation.  Similarly ,  ..  determines 
the  change  in  the  frequency  of  the  bifurcating  solution.  The  perturbational 
expansion  (15)  provides  a  good  approximation  to  the  periodic  solution  for 
computational  purposes. 

The  stability  of  the  bifurcating  oscillation  can  be  based  on  a  study  of  the 
Floquet  exponents  as  discussed  by  Iooss  and  Joseph  [34] .  The  essential  result 
is  that  the  periodic  solution  will  be  locally  stable  near  the  bifurcation  point 


if  the  eigenvalues  of  L  ,  other  than  ±iu>  ,  have  negative  real  parts  and  if 

U0  ° 

Re  j  (ij  )  is  positive. 

The  relevant  bifurcation  information  can  be  extracted  if  ;.0  ,  0  and 

a  (un)  are  computed.  The  algorithm  for  computing  each  of  these  is  presented 


next. 


An  Algorithm  for  the  Hopf  Bifurcation  Formulas . 

The  nonlinear  two-point  boundary  value  problem 

F(vU,u)  =  0  ,  G(vJ,u)  =  0  (19) 

must  first  be  solved  at  a  point  where  the  linearized  problem  (14)  has  a  pair 

of  purely  imaginary  eigenvalues  iioj^  •  Thus,  the  Hopf  point  is  located  by  finding 
a  root  of  Re  a(y).  This  is  accomplished  by  using  the  QZ  algorithm  to  com¬ 

pute  the  eigenvalues  at  each  point  along  the  steady  state  response  curves  and  then 
employing  the  bisection  or  secant  method  to  locate  . 

*  *  ;i0  i 

Let  L  denote  the  adjoint  differential  equation,  and  G  (v  ,u  i*)f  the 
u0  *  V  0 
adjoint  boundary  conditions.  The  eigenvectors  Cq  and  £q  are  then  computed 

from  the  eigenvalue  problems 
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(20) 


L,  ;0  =  H  '  Gv(V  'VV  =  °  1 

u 

=  V^olV-0-  (21) 

These  eigenvectors  are  then  normalized  by  requiring 

<  ?0'  «o  >  =  1  and  <  V  4  }  =  1  •  (22) 

Here  we  have  introduced  the  complex  inner  product 

<  v,«  >  =  /b  0  •  i  dx  (23) 

a 

where  >  *  $  denotes  the  dot  product  of  the  vector  functions  0  and  $■  . 

A  sequence  of  three  linear  nonhomogeneous  boundary  value  problems  must  next 
be  solved.  These  are 

=  -F,(v  °'-0) 


\r  W1  =  "2Fw(v  '  ~o‘ W 

0 

=  ‘2Gw(V  °-0h0'V 


(\  ’  2iu01,W2  “  "FW(V  'VW 


Gv(v  °'uo'W2)  "  ~Gvv(V  °'yoiS'V  ‘ 

The  derivatives  F^v  and  appearing  in  the  above  problems  are  computed  by 


the  rule 

,2 

F  (vP,uj$,v)  =  ■  rMr-  (vU  +  i  <j,  +  6  ,  .  (27) 

vv  -jz^'jz  1  2  '<5  =o  =0 

With  these  computations  complete,  we  now  compute  a^(yQ),  and  from 

y 

=  (Fvv(V  °'*Joh0»  %->  +  FVV1<V  °^0^0)  '^0  >  (28) 


+  u,o.  (un)}  =  -  <F_(v  ,u0k0^0,  ;0),;0  > 


-<Fw(v  'VW'S  >  {29) 

-  <pw(v  0  .%\lQ,*2).CQ  >  . 

The  first  author  evaluated  the  above  expressions  by  using  the  Stepleman  finite 

difference  scheme  to  solve  equations  (24-26)  and  using  Simpson’s  rule  to  numerically 

integrate  the  inner  products.  The  second  author  solved  the  boundary  value  problems 
with  the  de  Boor  and  Weiss'  collocation  code  and  integrated  the  inner  products  via 

Gauss-Lobat to  quadrature.  Both  methods  worked  well  and  yielded  quite  comparable 

results  for  o  (u  )  ,  and  cj  .  (a  program  for  the  Hopf  computations  for  the 

general  problem  (6)  will  be  available  from  the  second  author.) 

We  note  that  for  computations  where  numerical  sensitivity  was  observed,  the 

sensitivity -could  be  eliminated  by  more  accurate  computations  of  the  steady  state 

U0  * 

solution  v  ,  the  eigenvectors  and  ,  and  the  parameter  value  y^  at 

which  Re  a(y^)  =  0  .  This  sensitivity  is  not  removed  by  increasing  the  accuracy 
of  the  yQ  calculation  for  a  specific  discretization  as  in  the  case  for  ordinary 

differential  equations.  Some  obvious  tests  for  accuracy  are  an  orthogonality 

*  —  * 

relation  for  and  £  ,  <  )  =  0  ;  the  size  of  Re  a(uQ);  and  perhaps 

more  importantly,  comparisons  of  the  eigenvalues  of  the  linearized  problem  (14) 
and  its  adjoint. 


RESPONSE  CURVES 


In  this  section,  we  present  some  of  the  more  interesting  results  obtained 
by  using  the  above  numerical  techniques  to  study  the  dynamics  of  the  tubular 
reactor  which  arise  from  varying  the  Damkohler  number.  It  is  to  be  emphasized 
that  this  is  not  a  complete  cataloging  of  the  dynamic  capabilities  of  the  reactor 
but  rather  a  presentation  of  the  more  illustrative  cases  found  in  our  investiga¬ 
tions,  Many  of  the  response  curves  are  not  found  in  the  CSTR  studies  and  are 
thus  somewhat  unexpected.  The  cases  that  are  presented  can  be  viewed,  in  our 
opinion,  as  giving  a  much  more  complete  picture  of  the  interplay  between  the 
instabilities  in  the  steady  states  and  the  oscillatory  dynamics  that  can  arise 
during  operation  of  the  reactor. 

For  a  particular  combination  of  the  system  parameters,  we  summarize  the 
multiplicity  and  stability  of  the  steady  states  and  oscillatory  states  on  a 
response  curve  with  the  Damkohler  number  as  the  abscissa  and  the  maximum  temper¬ 
ature  of  the  steady  state  or  periodic  solution  as  the  ordinate.  The  steady 
state  solution  branches  are  computed  by  using  Keller* s  modification  of  the  Euler- 
Newton  procedure;  the  stability  of  these  branches  and  the  location  of  the  bi¬ 
furcation  points  are  determined  by  an  eigenvalue  analysis.  The  Hopf  bifurcation 
formulas  (16-18)  are  then  computed  to  give  the  direction  of  bifurcation  and  the 
stability  of  the  bifurcating  periodic  solutions.  By  using  the  asymptotic 
formula  (15)  for  the  solution  near  the  bifurcation  point,  we  are  able 
to  start  tracing  the  stable  oscillations  as  the  Damkohler  number  is  varied  by 
solving  the  full  parabolic  partial  differential  equations  with  PDECOL,  a 
general  code  based  on  the  method  of  lines  and  collocation  using  B-splines  [34]. 

The  first  example  is  illustrated  in  Figure  1  which  corresponds  to  the 

parameter  values  p  =  P  =  5 ,  B  =  0.5,  y  =  25,  8  =  3.5,  and  0_  =  1  . 

e,  e  0 

h  m 

"or  all  values  of  the  Damkohler  number,  D  ,  the  steady  state  is  unique  with 
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DAMKOHLER  NUMBER 


exchanges  m  stabi lity  at  D  =  0.262  and  D  =  0.295.  At  the  value  D  =  0.295, 


stable  periodic  solutions  bifurcate  to  the  left/  while  at  D  =  0.262,  the  solu¬ 
tions  also  bifurcate  to  the  left  but  are  not  stable.  PDECOL  was  used  to  trace 
the  stable  oscillation  from  D  =  0.295  down  to  D  -  0.260  where  the  stable 
oscillations  cease  to  exist  and  the  time-dependent  solutions  converge  to  the 
stable  steady  state  directly  below.  (We  conjecture  that  the  stable  branch  of 
periodic  solutions  connect  with  the  unstable  branch  emanating  from  D  =  0.262.) 

The  response  curve  dynamics  associated  with  varying  the  Damkohler  number 
can  now  be  explained  for  the  case  depicted  in  Figure  1  .  For  D  close  to  zero, 
the  reactor  operates  in  stable  state  of  low  temperature  and  low  conversion  of 
reactant  A  .  As  D  is  increased,  the  steady  state  remains  stable  and  the 
steady  state  temperature  rises  slowly  until  D  passes  through  D  -  0.262.  A 
jump  occurs  at  this  point  in  the  temperature  and  concentration  profiles  into 


the  stable  oscillation  directly  above.  These  oscillatory  solutions  remain  until 
D  reaches  D  =  0.295  after  which  the  reactor  operates  in  a  stable  steady  state. 


To  lower  the  conversion  and  temperature  in  the  reactor,  the  Damkbhler  number  is 
now  decreased.  At  D  =  0.295  a  small  stable  oscillation  in  the  temperature  and 
concentration  profiles  begins  to  appear.  These  oscillations  continue  until 
D  =  0.260  at  which  point  the  oscillations  disappear  through  a  jump  back  down 
to  the  lower,  stable  steady  state.  Both  of  these  jumps  may  be  thought  of  as 
ignition  or  extinction  processes,  respectively. 

In  Figure  2,  r,  is  decreased  to  2.5  and  a  region  of  three  steady  states 
appears.  Exchanges  of  stability  again  occur  at  the  Hopf  points,  and  at  the  upper 
point  (Li  =  0.1818),  the  behavior  of  the  oscillation  is  similar  to  that  of  Figure  1. 
The  stable  periodic  solution  bifurcates  to  the  left  and  its  amplitude  quickly  in¬ 
creases.  The  orbit  at  the  lower  Hopf  point  (D  =  0.165)  bifurcates  sharply  upward 
to  the  right  and  is  stable.  The  continuation  of  these  stable  branches  is  quite 
interesting.  Attempts  to  continue  the  larger  amplitude  periodic  branch  below 


D  =  0.172  failed  with  the  time-dependent  calculations  converging  to  the  lower 
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States  (Pe=5,  B=0.50,  y=25 


amplitude  periodic  solutions  regardless  of  the  initial  Conditions.  When  the  un¬ 
steady-state  calculations  were  started  on  the  lower  portion  of  the  branch  and  ; 
was  increased  above  0.172,  the  periodic  orbit  jumped  to  the  upper  solutions.  Tht 
temperature  profiles  from  the  upper  and  lower  periodic  branches  are  presented  in 
Figures  3  and  4.  The  explanation  of  this  seemingly  new  jump  phenomena  along  the 
periodic  branch  is  the  subject  of  current  investigation .  Since  the  response- 
curve  in  this  case  is  similiar  to  the  curve  in  Figure  5  ,  we  forego  a  discus  si.'  n 
of  the  possible  reactor  dynamics  for  this  case. 

The  dimensionless  heat  transfer  coefficient  is  lowered  to  2.35  in  Figure  r 
where  1-3-1-3-1  multiplicity  arises.  If  D  starts  near  zero  and  is  increased, 
the  reactor  temperature  increases  and  remains  in  a  stable  steady  state  until  the 
lower  limit  point  is  encountered  at  which  point  the  temperature  jumps  into  the 
much  higher  stable  oscillation.  This  stable  oscillation  surrounds  the  unstable 
steady  states  and  grows  steadily  until  a  jump  to  an  even  higher  stable  oscilla¬ 
tion  occurs  at  about  D  =  0.159.  The  amplitude  of  the  temperature  oscillation 
continues  to  grow,  peaks  out,  and  then  rapidly  decays  to  the  steady  state  at 
D  =  0.166.  The  reactor  can  now  be  extinguished  by  decreasing  the  Damkbhler 
number.  The  amplitude  of  the  periodic  orbits  bifurcating  at  D  =  0.166  grows 
rapidly,  peaks  out,  decays,  jumps  down,  continues  to  decay,  and  then  disappears 
as  the  time-dependent  solutions  converge  to  the  lower,  stable  steady  states 
beginning  at  the  lower  turning  point.  (It  should  be  pointed  for  this  case  that 
the  frequency  of  the  oscillation  decreases  along  the  lower  periodic  branch  as 
D  is  decreased.  The  period  of  the  oscillatory  solution  above  the  lower  limit 
point  is  approximately  five  times  greater  than  the  period  of  the  orbit  of  the 
Hopf  point,  D  =  0.166.) 
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igure  5.  Existence  of  Stable  Oscillatory  Solutions  Surrounding  Multiple  and  Unique 
Steady  States  (Pe=5,  B=0.50,  y=25,  and  6=2.35) 


The  existence  of  five  steady  states  for  the  reactor  is  shown  in  Figure  6  where 
is  decreased  to  2  »  A  Hopf  point  is  located  along  the  upper  section  of  the 
steady  state  branch,  and  an  unstable,  periodic  orbit  bifurcates  to  the  right  of 
this  point.  While  it  is  conceivable  that  a  turning  point  could  join  this  unstable 
branch  with  a  stable  one,  time-dependent  calculations  did  not  locate  any  stable 
period  solutions  in  this  example.  Furthermore,  u ,  (which  governs  the  asymptotic 
change  in  frequency)  is  a  large  negative  number  and  indicates  the  period  of  the 
unstable  solution  becomes  infinite  near  the  Hopf  point. 

It  is  a  straightforward  procedure  to  trace  the  bifurcation  diagrams  for  £  =  1 
and  3=0  to  complete  this  series  of  examples.  However,  Hopf  bifurcation  does 
not  appear  for  these  cases.  For  6=1,  the  reactor  exhibits  five  steady  states 
with  an  exchange  of  stability  at  each  limit  point  and  the  well-known  adiabatic 
result  of  three  steady  states  is  found  when  6=0. 

In  the  last  two  figures,  we  show  periodic  solutions  bifurcating  from  new 
patterns  of  multiple  steady  states  [24].  Five  steady  states  are  illustrated  in 
Figure  7  (pe  =  1,  B  =  0.50,  y  =  75,  and  6=4),  and  we  note  the  existence  of 
three  Hopf  points.  The  periodic  solutions  bifurcating  to  the  right  and  to  the 
left  from  the  intermediate  steady  states  are  unstable,  because  an  exchange  of 
steady  state  stability  does  not  accompany  the  bifurcation.  The  solutions  at  the 
upper  Hopf  point  bifurcate  to  the  left,  are  stable,  but  become  infinitely  periodic 
very  near  the  bifurcation  point. 

The  reactor  model  exhibits  one,  three,  five,  and  seven  steady  states  as 
well  as  bifurcating  periodic  solutions  in  Figure  8  (pe=l,  B=0.50,  y=125,  and 
3=4).  Both  of  the  periodic  orbits  bifurcate  to  the  left  and  both  become  infinitely 
periodic  near  the  Hopf  point.  Only  the  periodic  solutions  bifurcating  just  to 
the  right  of  the  upper  quench  point  are  stable. 
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igure  8.  Existence  of  Multiple  Periodic  Solutions  and  Seven  Steady  States  (Pe— 1,  B-0.50,  y-125r  and 


In  Figures  6-8,  the  response  curves  become  quite  complex  exhibiting 
several  steady  states  and  oscillatory  solutions.  However,  the  dynamic  capa¬ 
bilities  of  the  reactor  in  these  examples  are  somewhat  limited  since  many  of 
the  states  are  unstable  over  a  large  range  of  the  Damkbhler  number.  Ingition 
is  characterized  by  increasing  D  past  a  lower  limit  point  so  that  the  reactor 
jumps  from  a  stable,  low  temperature  steady  state  to  a  stable,  high  temperature 
steady  state.  The  reactor  is  extinguished  by  decreasing  D  until  a  jump 
occurs  from  an  upper  stable  state  (either  a  steady  state  or  a  periodic  state 
with  small  amplitude)  down  to  a  lower,  stable  steady  state.  In  these  last  three 
examples,  jumps  from  lower  steady  states  to  large  amplitude  periodic  solutions, 
jumps  between  oscillatory  states,  and  possible  periodic  reactor  operation  over 
large  regions  of  the  Damkohler  number  are  not  possible. 
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CONCLUSIONS 


We  have  presented  numerical  bifurcation  techniques  which  can  determine 
all  possible  steady  and  stable  oscillatory  solutions  exhibited  by  distributed 
parameter  models.  These  techniques  were  applied  to  the  nonadiabatic  tubular 
reactor  and  illustrated  several  types  of  steady  state  and  oscillatory 
phenomena.  Many  of  the  results  were  anticipated  from  studies  of  the  CS TR, 
but  it  is  clear  that  all  the  possible  phenomena  exhibited  by  the  tubular 
reactor  has  not  been  uncovered.  The  effects  of  parameters  such  as  the  feed 
temperature  and  flow  velocity  which  could  be  used  to  control  the  reactor  are 
largely  unknown.  With  the  proper  motivation  and  support,  our  techniques  can 
be  used  to  solve  these  problems.  Furthermore,  because  of  their  generality, 

' 

these  numerical  methods  can  also  be  applied  to  a  broad  array  of  models  found 
in  combustion  theory  and  mathematical  biology  as  well  as  in  chemical  reactor 
theory. 


i 
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NOTATION 


A 

B 


D 

D 

e 

E 

AH 

k 

e 

L 

P 

Pe 

n 

peh 

R 

t 

T 

T 


cross-sectional  area  of  the  reactor,  m 
frequency  factor,  s 

dimensionless  heat  of  reaction,  A  H  c_  /pC  T 

u  p  0 

concentration,  mol/m^ 

3 

inlet  concentration ,  mol/m 
specific  heat,  J/mol°K 

— v 

Damkohler  number,  A  e  L/v 

2 

effective  diffusivity,  m  /s 

activation  energy,  J/mol 

heat  of  reaction,  J/mol 

effective  thermal  conductivity  j/s  m°K 

reactor  length ,  m 

reactor  perimeter,  m 

Peclet  number  for  mass  transfer  vL/D^ 

Peclet  number  for  heat  transfer  p C  vL/k 

P  e 

universal  gas  constant 
time ,  s 
temperature,  °K 
inlet  temperature,  °K 


dimensionless  axial  distance,  x/L 

2 

heat  transfer  coefficient,  j/m 
velocity,  m/s 
axial  distance,  m 
dimensionless  concentration ,  c/c 


s  °K 
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